#include "math.h"
#include <algorithm>
#include <map>
#include <string>
#include <sstream>
#include <fstream>
#include <iostream>
#include <vector>
#include <map>
#include <sys/time.h>
#include <unistd.h>
#include <limits>
#include "bessel.cpp"


using namespace std;

#define PERIOD 3
const char ALPHABET[]={'a', 'c', 't', 'g'};
#define ALPHABET_SIZE 4

typedef map<char, int> ICol;
typedef vector<ICol> IMatrix;
typedef map<char, float> FCol;
typedef vector<FCol> FMatrix;

typedef pair<float, int> DPCell;
typedef vector<DPCell> DPCol;
typedef vector<DPCol> DPMatrix;

struct DPResult {
    float reward;
    int right_x;
    int right_y;
};

struct DPResultExt {
    float reward;
    int right_x;
    int right_y;
    int left_x;
    int sh_num;
    int sh1;
};

struct CountShiftsRes {
    int sh_num;     // число сдвигов фазы
    int sh1;        // координата первого слева сдвига
    int left_x;     // координата левого конца найденного участка
};

template<typename ItemType>
ItemType accumulate(vector<ItemType> vec) {
    ItemType sum=0;
    for (int i=0; i<vec.size(); i++) {
        sum += vec[i];
    }
    return sum;
};

extern "C"
void print_imx(IMatrix mx) {
    int i,j;
    for (j=0; j<ALPHABET_SIZE; j++) {
        for (i=0; i<PERIOD; i++)
            printf("%5i\t", mx[i][ALPHABET[j]]);
        printf("\n");
    }
};

extern "C"
void print_fmx(FMatrix mx) {
    int i,j;
    for (j=0; j<ALPHABET_SIZE; j++) {
        for (i=0; i<PERIOD; i++)
            printf("%.3f\t", mx[i][ALPHABET[j]]);
        printf("\n");
    }
};

extern "C"
IMatrix create_mx(string seq) {
    int i,j;
    char c;
    IMatrix mx;

    // init mx by 0
    for (i=0; i<PERIOD; i++) {
        ICol col;
        for (j=0; j<ALPHABET_SIZE; j++) {
            col[ALPHABET[j]] = 0;
        }
        mx.push_back(col);
    };
    
    // fill mx by seq
    for (int i=0; i<seq.length(); i++) {
        mx[i%PERIOD][seq[i]]++;
    };
    return mx;
};

extern "C"
FMatrix normalize_mx(IMatrix mx) {
    int i,j,m;
    float p;
    
    vector<int> y;
    for (i=0; i<PERIOD; i++) {
        m = 0;
        for (j=0; j<ALPHABET_SIZE; j++) {
            m += mx[i][ALPHABET[j]];
        };
        y.push_back(m);
    }
    
    vector<int> x;
    for (j=0; j<ALPHABET_SIZE; j++) {
        m = 0;
        for (i=0; i<PERIOD; i++) {
            m += mx[i][ALPHABET[j]];
        };
        x.push_back(m);
    }
    
    float L = float(accumulate(y));
    
    FMatrix nmx;
    for (i=0; i<PERIOD; i++) {
        FCol col;
        for (j=0; j<ALPHABET_SIZE; j++) {
            p = y[i] * x[j] / L / L;
            //printf("NMX: %i %f %f\n", mx[i][ALPHABET[j]], p*L, sqrt(L * p * (1.0 - p)));
            if (p > 0)
                col[ALPHABET[j]] = (mx[i][ALPHABET[j]] - p * L) / sqrt(L * p * (1.0 - p));
            else 
                col[ALPHABET[j]] = 0.0;
        }
        nmx.push_back(col);
    };
    //print_imx(mx);
    //print_fmx(nmx);
    return nmx;    
};

extern "C"
FMatrix create_norm_mx(string seq) {
    return normalize_mx(create_mx(seq));
};

extern "C"
float Bessel2Norm(float x) {
    int ai = abs(int(x*10));
    if (ai > 1000)
        ai = 1000;
    //printf("%f\t%f\n", x, BesselNormValues[ai]);
    if (x>=0)
        return BesselNormValues[ai];
    else
        return -BesselNormValues[ai];
};

extern "C"
float norm_dist_mx(FMatrix nmx1, FMatrix nmx2) {
    int i,j,x;
    float D = 0.0;
    for (i=0; i<PERIOD; i++) {
        for (j=0; j<ALPHABET_SIZE; j++) {
            D += Bessel2Norm( nmx1[i][ALPHABET[j]] * nmx2[i][ALPHABET[j]] );
        };
    };
    return D;
};

extern "C"
CountShiftsRes count_shifts(DPMatrix &dp_mx, int max_i, int max_j, bool local=true) {
    int shift = max_j;
    int cnt_shifts = 0;
    int sh1 = 0;
    int i;
    
    //printf("max_i = %i; max_j = %i\n", max_i, max_j);
    for (i=max_i; i>=0; i--) {
        DPCell& cell = dp_mx[i][shift];
        if (local and cell.first <= 0)
            break;
        //printf("shift = %i\n", shift);
        if (cell.second != shift) {
            sh1 = 0;
            cnt_shifts ++;
        }
        shift = cell.second;
        sh1 ++;
    }    
    return {cnt_shifts, sh1-1, i+1};
}

extern "C"
DPResult fill_dp_mx(DPMatrix &dp_mx, int reward, int penalty, int win_length, string base_seq, string seq, bool local=true) {
    int i, j, k, x, _max_j;
    int max_x = 0;
    int max_j = 0;
    float w, max_w;
    float max_reward = 0;
    int length = seq.length();
    float MIN_W;
    
    if (local)
        MIN_W = 0.0;
    else
        MIN_W = numeric_limits<int>::min();

    IMatrix base_imx = create_mx(base_seq);
    //print_imx(base_imx);
    FMatrix base_mx = normalize_mx(base_imx);
    //print_fmx(base_mx);

    vector<float> deltas;
    deltas.resize(PERIOD);
    
    DPCell max_cell;
    
    int stat_0 = 0;
    int stat_all = 0;
    
    vector<FMatrix> bmxs;
    
    for (i=0; i<PERIOD; i++) {
        FMatrix temp_mx;
        for (j=0; j<PERIOD; j++) {
            temp_mx.push_back(base_mx[(j+i) % PERIOD]);
        };
        bmxs.push_back(temp_mx);
    };

    DPCol col;
    for (i=0; i<PERIOD; i++) {
        col.push_back(make_pair(0, i));
    }
    dp_mx.push_back(col);
    
    for (x=0; x<length-win_length-1; x++) {
        //printf("%i ################################\n", x);
        //printf("deltas:  ");
        for (i=0; i<PERIOD; i++) {
            deltas[i] = norm_dist_mx(
                bmxs[x % PERIOD],
                create_norm_mx(seq.substr(x+i, win_length))
            );
            //printf(" %f ", deltas[i]);
        };
        //printf("\n");
        //print_fmx(bmxs[x%PERIOD]);
        //printf("***\n");
        //print_fmx(create_norm_mx(seq.substr(x, x+win_length)));
        //printf("***\n");
        
        
        DPCol col;
        for (i=0; i<PERIOD; i++) {
            max_w = 0.0;
            _max_j = i % PERIOD;
            for (k=i; k<i+PERIOD; k++) {
                j = k % PERIOD;
                if (i==j)
                    w = dp_mx[x][j].first + deltas[j] - reward;
                else
                    w = dp_mx[x][j].first - penalty;
                if (w > max_w) {
                    max_w = w;
                    _max_j = j;
                };
            };
            col.push_back(make_pair(max_w, _max_j));
            if (max_w > max_reward) {
                max_reward = max_w;
                max_j = _max_j;
                max_x = x;
            }
        }
        dp_mx.push_back(col);
        //printf("%5i: ", x);
        //for (i=0; i<PERIOD; i++) {
            //printf("(%7.2f, %i), ", col[i].first, col[i].second);
            //stat_all ++;
            //if (col[i].first == 0)
                //stat_0 ++;
        //}
        //printf("\n");
    };
    //printf("%i / %i = %f\n", stat_0, stat_all, float(stat_0)/stat_all);
    
    if (local) {
        return {max_reward, max_x, max_j};
    }
    else {
        max_x = dp_mx.size()-1;
        max_j = -1;
        max_reward = MIN_W;
        DPCol last_col = dp_mx[max_x];
        for (i=0; i<PERIOD; i++) {
            if (last_col[i].first > max_reward) {
                max_reward = last_col[i].first;
                max_j = last_col[i].second;
            }
        }
        return {max_reward, max_x, max_j};
    }
};


extern "C"
DPResultExt create_dp_mx_c(int reward, int penalty, int win_length, char *base_seq_c, char *seq_c, bool local=true) {
    string s1 = (const char*) base_seq_c;
    string s2 = (const char*) seq_c;
    DPMatrix dp_mx;
    
    DPResult fres = fill_dp_mx(dp_mx, reward, penalty, win_length, s1, s2, local);
    CountShiftsRes cres = count_shifts(dp_mx, fres.right_x, fres.right_y, local);
    return {
        fres.reward,
        fres.right_x,
        fres.right_y,
        cres.left_x,
        cres.sh_num,
        cres.sh1,
    };
}


